FAMILIES OF SURFACE GAP SOLITONS AND THEIR STABILITY VIA THE 
NUMERICAL EVANS FUNCTION METHOD 
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Abstract. The nonlinear Schrodinger/Gross-Pitaevskii equation with a linear periodic potential and a nonlincarity 
coefficient V with a discontinuity supports stationary localized solitary waves with frequencies inside spectral gaps, so 
called surface gap solitons (SGSs). We compute families of ID SGSs using the arclcngth continuation method for a 
range of values of the jump in T. Using asymptotics, we show that when the frequency parameter converges to the 
bifurcation gap edge, the size of the allowed jump in T converges to for SGSs centered at any x c € K. 

Linear stability of SGSs is next determined via the numerical Evans function method, in which the stable and 
unstable manifolds corresponding to the solution of the linearized spectral ODE problem are evolved up to a common 
location where the determinant of their bases, i.e., the Evans function, is evaluated. Zeros of the Evans function coincide 
with eigenvalues of the linearized operator. Far from the SGS location the manifolds arc spanned by exponentially 
decaying/increasing Bloch functions. As we show, evolution of the manifolds suffers from stiffness. A numerically 
stable formulation is possible in the exterior algebra formulation and with the use of Grassmanian preserving ODE 
integrators. Eigenvalues with a positive real part larger than a small constant arc then detected via the use of the 
complex argument principle and a contour parallel to the imaginary axis. The location of real eigenvalues is found 
via a straightforward evaluation of the Evans function along the real axis and several complex eigenvalues are located 
using Miillcr's method. The numerical Evans function method is described in detail in order to facilitate its use as 
a practical tool for locating eigenvalues. Our results show the existence of both unstable and stable SGSs (possibly 
with a weak instability), where stability is obtained even for some SGSs centered in the domain half with the less 
focusing nonlincarity. Direct simulations of the PDE for selected SGS examples confirm the results of Evans function 
computations. 
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function method, exterior algebra, argument principle, Miiller's method 
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1. Introduction. Waves in periodic media with a nonlinear response are of fundamental im- 
portance in several areas of physics including propagation of light in nonlinear photonic crystals and 
evolution of Bose Einstein condensates (BECs) loaded on optical lattices. A classical and highly 
universal one dimensional model in these settings is the periodic nonlinear Schrodinger (PNLS) / 
Gross-Pitaevskii equation 

id t ip(x,t) +dli>{x,t) - V(x)ip(x,t) +T(x)\i){xit)\ 2 iP(x,t) = (1.1) 

on x G R, t > with V(x + d) = V(x) V x g R. Equation f) 1 . 1 [) is a model for BECs confined via a 
magnetic or optical trap (in the directions transversal to x) to a quasi-lD geometry and loaded on an 
optical lattice [9l [15] and also for optical beams propagating in Kerr nonlinear media in a direction 
orthogonal to the x— axis along which the linear part of the refractive index of the medium is periodic 
[t7] . In the latter case t is a spatial variable along the propagation direction. We investigate stationary 
solitary waves for nonlinear interfaces 

r(as) = r+ x[o,oo)(x) + r_ X(-oc,o)( x )> where r ± e R - ( L2 ) 

In optics this coefficient describes two media with different values of the cubic susceptibility affixed 
at x = while in BECs it corresponds to two condensates with different s-wavc scattering lengths. If 
T(x) > 0, the BEC interaction is 'attractive' at x; in optics the medium is called 'focusing' at x. In 
the opposite case the respective terms are 'repulsive' and 'defocusing'. 
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Localized solitary waves are of particular physical as well as mathematical interest due to their 
properties of finite energy and constant shape throughout evolution. In optics they are considered as 
viable candidates for bit carriers. Solitary wave solutions ip(x,t) = e~ ltlt cj)(x), /i € K of with a 
periodic V and with exponentially localized <f> are called surface gap solitons (SGSs) when F + =/= F_ in 
(|1.2[) and gap solitons (GSs) when T is periodic. The spatial profile <fi satisfies the stationary periodic 
nonlinear Schrddinger equation 



The word 'gap' in the names of GS and SGS comes from the fact that the frequency parameter /i 
must lie in a gap of the spectrum o~(L) of the operator 



(S)GSs arise due to a balance between the linear periodicity induced dispersion and the nonlinear 
focusing or defocusing. GSs of (jl.ip were studied theoretically, for instance, in [171131]. In [ISI 
existence of gap soliton ground states in all gaps of cr(L) for the d-dimensional PNLS was proved 
rigorously using variational methods. GSs have been also observed experimentally in optics [7J [TB] as 
well as in BECs [16]. 

One dimensional SGSs have been previously studied mainly at linear interfaces, i.e., for T = const, 
and V(x) = V\{x) X[o,oo){x) + V-z{x) X(-oo,o)( x ): ^1,2(2: + d\ , 2 ) = Vi. 2 {x), see the numerical studies 
e.g. in [35] [Mj. Analogous linear interfaces have been considered in 2D J2S][27]. SGSs have been 
observed also experimentally in photonic crystals [331 HH] • 

At the nonlinear interface p. 21) SGSs were studied via numerical and asymptotic methods in [13] . 
Using the implicit function theorem, their existence was shown to hold via a parameter continuation 
from GSs with L =const. for sufficiently small \T + — L_| under a non-degeneracy condition on the 
GS. Families of SGSs were then numerically continued in F_ for L + fixed up to a critical point where 
the corresponding Jacobian becomes singular. Here we first use the numerical arclength continuation 
to extend these families past the critical points, which are all apparently simple folds of the SGS 
families. The families are continued up to a chosen threshold value of T_ or total energy. Next, we 
investigate linear stability of these SGSs via the numerical Evans function method and identify stable 
and unstable segments of the SGS families. We show that folds are often associated with a stability 
change. 

The phenomenon of localization in a medium with a purely nonlinear interface relies crucially 
on the presence of a linear structure. Without a linear structure, i.e., for V = const., no localized 
solutions of (|1.3|) exist as shown in the following lemma. This result also appears (for real <j)) in Sec. 



Lemma 1.1. If V = const, and L + 7^ F_, then the stationary periodic nonlinear Schrddinger 
equation Q1.3P has no nontrivial solutions in 

Proof. Due to the presence of the term [i§ in (|1.3[) we can assume without any loss of generality 
that V = 0. As a first order system in real variables (/)r c := Re(0) and <fii m := lm(0) equation (|1.3p 
reads 



which is a Hamiltonian system with H (</> Rc , Im , <^ c , $ m ) = \ [(<t>' Rc ) 2 + (</>{ m ) 2 + m(0r c + ^D] + 

jF(x) ((/)^ c + 0j m ) . Because F(x) = F + X[o,oc)( x ) + T- X(- O o,o)( a; ); the Hamiltonian H is conserved 
on each x > and x < 0. As e iJ 1 (R), both 4>(x) and (f>'(x) decay to as \x\ — > 00 so that 
H (4>n e (x), 0i m (a;), 0^ e (x), <p[ m (x)) -> as |a;| -> 00, and thus H = 0. 



<f)"(x) + n<f>{x) - V{x)(j){x) + T(x)\cj)(x)\ 2 (j)(x) =0, 



(1.3) 



L:= -d 2 x + V{x), iei 



3.3 of [13]. 
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Differentiation yields 



o = -^h = <j)' Kc [< e + /i0 Rc + r( x ) (0L + 0L) 0R C ] + 

1 2 
01m Km + A^Im + r(x) ((f) Im + Rc ) 0I m ] + -(!+- r_)6(x) ((£ Re + <^ m ) , 



(1.4) 



where 5(:r) is the Dirac delta. Integrating (fl~4]) over E, we get = j(T + - r_) (0 Re (O) + 0? ro (O)) so 
that ^ R e(0) = <t>i m (0) = 0. 

Finally, due to the C 1 -nature of <f> we have also = lim x ^o+ ^(0Re(^), 0im(^), 0R e ( x ); 01m ( x )) = 



5 ((0Re(O)) 2 + (0L(O)) 2 ) • In summary Rc (O) = = 0i m (O) = ^ m (0) - 0. Hence = 0. □ 



As a result, solitons of the constant coefficient nonlinear Schrodingcr equation cannot be continued 
in (T_ — r+) from T_ = T+. This is due to the presence of the spatial shift invariance, which renders 
the corresponding Jacobian singular so that a parameter continuation fails. A periodic linear structure 
can thus be viewed as one of the simplest structures supporting localized waves at nonlinear interfaces. 
In the rest of the paper we limit our attention to real SGSs profiles cf>(x) : E — > M. These satisfy 



so that the first variation of (</>) is equivalent to (|1.5[) . 

Understanding stability of SGSs with respect to perturbations of the initial data 4>(x) is of crucial 
importance for predicting physically and numerically observable solutions since experiments as well 
as numerical simulations generate unavoidable error. We investigate linear stability by inspecting the 
spectrum of the corresponding linearized operator. Standard methods for numerically approximating 
spectra of operators include direct eigenvalue computations in a finite dimensional approximation, e.g. 
finite differences or finite elements; Floquct-Bloch theory, where the problem is artificially treated 
as periodic [TT]; and the numerical Evans function method [12j |4j HTJ [19], out of which the last 
method is in principle restricted to one dimensional problems. Direct eigenvalue computations in a 
finite dimensional space typically suffer from spurious eigenvalues, which are difficult to identify |12j . 
Although certain techniques for tackling this spectral pollution are available [KB [35] , we choose the 
Evans function method, which does not suffer from spectral pollution. A further advantage of the 
Evans function method is the possibility to determine the number of eigenvalues within a region of the 
complex plane with a single contour integration via the use of the argument principle. Nevertheless, 
we do not claim the superiority of this method over others. In the Evans function approach finding 
eigenvalues is reformulated as determining linear dependence of the stable and unstable manifolds of 
the zero solution of the linearized ODE. In the end this problem reduces to the evolution of an ODE 
system and the computation of a determinant, called the Evans function. Zeros of the Evans function 
coincide with eigenvalues of the linearized operator. As we show, a numerically stable and accurate 
evaluation of the Evans function, however, requires a careful computation of exponentially exploding 
solutions, the use of exterior algebra, which avoids numerical stiffness while preserving analyticity of 
the Evans function as well as the use of an ODE integrator which approximately preserves the weak 
Grassmanian invariant of the ODE system. Analyticity is essential in our method as we apply the 
argument principle in order to determine eigenvalues in the right half complex plane. 

The aim of the paper is twofold. Firstly, it is to investigate properties of SGSs of (jl.ip . including 
their stability, and secondly to develop a practical numerical Evans function method combined with 
the argument principle for determining linear stability in problems with periodic coefficients. At the 
same time we point out possible difficulties of the method. 



4>"{x) + n4>{x) - V{x)(j){x) + r(x)(j) 3 (x) = o, x e E. 

The C 1 -solutions of (|1.5[) are critical points of the total energy 



(1.5) 




(1.6) 
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The rest of the paper is structured as follows. In the remaining part of the introduction we 
give a brief overview of Floquct theory and define some quantities needed in later sections. Section 
[2] discusses the numerical method of arclength continuation as applied to computing SGS families. 
Results of these computations are then presented and discussed for SGSs in the first two spectral 
gaps. Section [3] is devoted to linear stability of SGSs. After summarizing available information 
about the spectrum of the linearization we describe the Evans function method for locating the point 
spectrum. Next, we discuss difficulties in evaluating the Evans function in a numerically stable and 
accurate way and offer their solutions. Results of Evans function computations are then presented for 
the SGS families computed in Section [21 Finally, direct numerical simulations of the PDE (jl.lj) arc 
performed for selected SGS examples to demonstrate their (un)stable evolution. 

1.1. Review of Floquet Theory. The linear part of the stationary equation (|1.5I) is Hill's 
equation 

Lq = -dlq + V(x)q = fiq, x£ R. (1.7) 

An overview of Floquet theory for Hill's equation can be found in [33[ I14j. The spectrum <j(L) 
can be determined by studying the trace of the monodromy matrix of the first order version of (|1.7|) . 
The spectrum is composed of bands, a(L) = Ufc e N[s2fc-i, S2fc], where S2k+i > S2fc > S2fe-i- We define 
int(cr(L)) = Ufe e N(s2fc-i, S2fc) and da(L) = {si, s 2 , ■ ■ .}■ For fi £ C \ da(L) there are two linearly 
independent solutions of (|1.7p . so called Bloch waves, which have the form 

qx{x) = p!{xy hx , q 2 {x) = p 2 {x)e-' lkx , 

where pia{ x + d) = pi,2{x) Vs G R and k £ C. The periodic parts pi.2 satisfy 

- {d x ± ik) 2 Ph2 (x) + V(x) Ph2 {x) = W , 2 (x), x £ [0, d] (1.8) 

with periodic boundary conditions. 

For ix £ int(cr(i)) the Bloch waves are bounded since k £ R. For fi £ da(L) one solution is a 
periodic Bloch wave, and thus remains bounded, while the other linearly independent solution grows 
linearly in x. We denote the periodic Bloch function at the edge [i = s n £ dcr(L) by q n . When 
/i £ C \ da(L) 7 the Bloch waves grow exponentially in either x direction since Im(fc) =/= 0. 

2. SGS Existence and Numerical Continuation. In [13] an implicit function theorem argu- 
ment was used to show that gap solitons </>gSj i-e. solutions of (( 1 . 5|) with r = Tq, can be continued to 
SGSs as r_ departs from T + = To provided the Jacobian J = —d% + V(x) — 3Tq(/)q S (x) is nonsingular. 
This condition on J is expected to hold due to the absence of shift invariance in equation (| 1 . 5|) . Note 
that the continuation result holds for GSs centered at any x c £ R. In |13j continuation in T_ is 
carried out numerically for x c = and V(x) = sin 2 (7ra;/10) up to a critical value of T_, where J 
becomes singular. Here we show using the numerical arclength continuation method [301 129| that the 
SGS families simply turn back with respect to (r_ — T+) at this critical point and can be continued 
further. Often numerous turning points occur within an SGS family. These results suggest that the 
critical points are simple folds, where the kernel of J is only one dimensional. Arclength continuation 
has been used in the context of the discrete NLS for instance in [35] ■ 

Let us denote equation flL5J as G((f>, T_) = with T + fixed and (f> and T_ unknown. Due 
to the jump in the coefficient T the solution is only C 1 regular and we thus discretize (|1.5[) via 
H 1 conforming finite elements (using a FEM package by M. Richter from the Karlsruhe Institute of 
Technology). We use elements of 3rd order (p — 3) on a domain [-R, R] with R > large enough so 
that a given 0gs is well decayed at x = ±R (typically </>gs ~ 10~ 14 for x £ [-R, — jqR] U [^R, R})- 
We use the element size h < 0.06 in all computations. As we show in Sec. 13.5. It the accuracy of 
the Evans function, measured by how well the Evans function recovers the zero eigenvalue of the 
linearization, improves as h — > 0. 
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In order to describe the arclength continuation, let us denote the resulting FEM-discrctized system 
of N algebraic equations by 

G(0,r_) = O, (2.1) 

where </> = (4>k)k=i and <\>k is the approximation of <j>(xk),Xk = —R + {k — k e {1, . . . ,N}. In 

the arclength continuation method the solution family (</>,r_) is parametrized by the arclength r of 
the curve generated by the family in the (||0||,r_) plane with a selected norm || ■ ||, e.g. the ^-norm. 
System (|2.1j) is underdetermined and the arclength method appends it with an equation that ensures 
continuation in the tangent direction. Namely, to extend the family from r = To to r = t\ , one solves 

G(0,r_) = o 

(<?(n) - 0(r o )) T J^(r ) + (r_(n) - r_(r )) ^-(r ) - (n -r ) = 

for the vector ((jFiji), T_ (ri)) <E R JV+1 , which we do by Newton's iteration using ((^(tq), T _(t )) + (t 1 — 
r o)(^:^ r ( T o): ^r_(ro)) as the initial guess. After convergence the tangent vector at r = T\ is ob- 
tained by solving 



Gjin) Gr_(n) 
Mr ) fr_(r ) 



where G^ is the N x N Jacobian matrix ^J^-, ■ • ■ , -§^j and Gr_ = The solution vector is 

normalized so that || J^0(ti)||^ 2 + (£f^-( T i)) 2 = 1- The first continuation step from the GS at r = 0, 
where the tangent vector is not available, is performed via a standard continuation in the parameter 
r_ . Newton's method is terminated when the L 2 -norm of the residual decreases below a tolerance 
(typically 10" 10 ). 

In |41j it is shown that families of GSs originate at upper edges of spectral gaps for the focusing 
nonlincarity To > and at lower gap edges for To < 0. They bifurcate from the zero solution at the 
corresponding edge, which we call the bifurcation edge. For even potentials V with one maximum and 
one minimum in each period there are two distinct families of GSs, namely those symmetric about the 
minimum location, so called onsite GSs, and those symmetric about the maximum location, offsitc 
GSs [41j . GSs centered at x = can be computed via a numerical continuation in the parameter [i 
starting near a selected bifurcation gap edge fx = s n . As an initial guess for GSs near s„ we use the 
slowly varying envelope approximation eA n (X)q n (x), X = ex, where < e 2 << 1 is the distance 
of /i to the gap edge: \i = s n + fie 2 and fl = ±1 for lower and upper gap edges, respectively. q n is 
the periodic Bloch function at fj, = s n and A n satisfies the constant coefficient homogenized nonlinear 
Schrodinger equation, see e.g. j5], Sec. IV in [41] or Sec. 3.3 in [13] . 

OA* + vA' 7 [ + pTA 3 = 0, (2.3) 

where v = 1 + 2(q' n , q n )L 2 (o,d) and p = U^nH^o d ^ with q n being a (periodic) 'generalized Bloch func- 
tion'. We use the explicitly known ground state solution of (|2.3p . namely A n (X) = ^/^Pscch ^ =p-X 

In this sense we construct an approximation to 'fundamental' GSs. Choosing an excited state solution 
to (|2.3p . a different family of GSs would result. 

2.1. Numerical Continuation Results. We restrict our attention to the semi-infinite gap 
(— oo, si) and the first finite gap (s2, S3) and first compute families of both on- and offsitc GSs bifur- 
cating from upper gap edges si and S3 for T = 1 and from the lower edge S2 for T = — 1. In all of the 




Fig. 2.1. Families of onsite gap solitons in the semi-infinite and the first finite gap for the potential (1-'.4t (with 
xq = 0). (a) GS continuation curves in the frequency - total power coordinates; (b)-(d) GS profiles for the points 
labeled in (a). 



numerical computations we use 

V{x) = sin 2 U^^j , d = 10, (2.4) 

where Xq = for onsite and Xq = d/2 = 5 for offsite GSs. The shift Xq is introduced in order to ensure 
that both onsite and offsite GSs are centered at x = 0, the interface location for the subsequent SGS 
study. The first 10 edges of the spectrum a(L) for (|2.4[) have been computed using the trace of the 
monodromy matrix |14| to be 

(si, s 2 , sio) « (0.2832, 0.2905, 0.7467, 0.8433, 1.0568, 1.4069, 1.4505, 2.0989, 2.1022, 2.9806). 

Note that the edges S4, . . . , sio are needed below in Figs. 13.31 Figs. 12.1( a) and 12.2( a) show families of 
both onsite and offsite GSs 4>gs m the (/z, ||</>gs 11^,2 ) plane. In optics ||</>|| 2 2 is called the total power 
of <j>. For the labeled points A on — H on and A Q s — H s the GS profiles are also plotted. The ^-values 
of the selected GSs are listed below. 
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Families of SGSs bifurcating from these GSs have been computed next via the arclength continu- 
ation described above holding r + fixed at the corresponding GS value, i.e., either r + = 1 or r + = — 1. 
The points A on ~ H on and A Q g — H s have been selected as representative examples. This set includes 
points near bifurcation edges, where ||0||l2 — > as \i approaches the edge, near the opposite edges, 
and also deep inside the gaps. Figs. 12.31 and 12.41 show these SGS families in the (r_,£^(^>)) plane 
for the onsite and offsite case, respectively. The choice of E^((f)) over ||</>||^ 2 has been made due to 
a clearer identification of fold locations in the former case. The folds are shown in Sec. 13.51 to often 
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Fig. 2.2. Families of offsite gap solitons in the semi-infinite and the first finite gap for the potential 1 12. 41 1 (with 
xo = d/2 = 5). (a) GS continuation curves in the Jreguency - total power coordinates; (b)-(d) GS profiles for the 
points labeled in (a). 

correspond to stability changes of the SGS and in all the considered cases they are accompanied with 
a change in the number of eigenvalues in the right half complex plane. The black dots in Figs. I2.3l and 
12.41 correspond to the GSs from which continuation was started in both the T_ < T + and T_ > T + 
directions. In Fig. 12.51 we plot for illustration the profiles of several SGSs along the family C ff. 

2.1.1. Observations on Numerical Results. 

Occurrence of Translated Bifurcation GSs. As one can see in Fig. 12.31 fa), (c), and (f), families of 
SGSs may return to the point corresponding to the bifurcation GS, i.e. (r_, E u (</>)) = (r + , -E m (^>gs))- 
In fact, we show in Fig. 12.61 (using the families C on and F on as examples) that the solutions cor- 
responding to all these intersections of (r+, i5 M (</>Gs)) are GSs. Each return to (r + , ^(c^qs)) is, 
however, a GS centered at a different extremum of V. We expect that generically for a bifurcation 
from an onsite GS 0gs centered at x = each half of the SGS family, i.e., the one bifurcating in T_ 
to the left of T + and the one bifurcating to the right, the n— th return to (r+, £^(0gs)) is an onsite 
GS centered at x = nd or x = —nd. This has been checked to hold in all the computed examples with 
SGS families that are open curves in the (T_, E^(4>)) plane. Note that it is impossible to conclude 
based on the equalities T_ = T + and E^((j)) — £^(</>gs) that the solution <fi must be a GS but our 
numerics suggest this. Fig. 12.61 demonstrates this for the cases C on and F on , where returns with n = 1 
occur for both. 

For closed SGS family curves, like the case E Q g in Fig. 12.41 this rule is not expected to hold. 
Upon return of the family E g to (r_,£^(0)) = (r+, £^(<^>gs)) the solution is, indeed, identical to 
the bifurcation GS 4>gs centered at x = 0. As Fig. 12.71 shows, the other solution that lies on the line 
T_ = r + and belongs to the E Q g family is also an offsite GS centered at x = but different from </>gs- 
This is possible due to non-uniqueness of bound states of (|1.5p . Part (b) of Fig. 12.71 plots the center 



. xicr 



(a) r-1 



E. BLANK AND T. DOHNAL 
(b) r=l 




(e) r-1 



-2 r -1 



1 r -0-5 o 




6 r - 4 



Fig. 2.3. SGS families bifurcating from the onsite GSs A on — H on in Fig. \2.1\ Note that the curves were computed 
up to |r_| = 30 or \Efi{4>)\ = 100, whichever occurred first, and the behavior is qualitatively the same as in the displayed 
windows. 



of "mass" 

xcm{4>) ■■= \ x\4>(x)\dx / / \4>{x)\dx 

JR I JR 

for the family and shows that the whole family E q r consists of solutions whose center of mass does 
not move far from the minimum point x = of V compared, e.g., to the cases C on , F a g below. 

It is an open question how many returns to (r_,£^(0)) = -Em(</>gs)) occur in a given SGS 
family. Note that none of the open families of SGSs in Fig. 12.41 which bifurcate from offsite GSs, 
returns to (r + , £^(0gs)) in contrast with the onsite case where these returns occur, e.g. for A on , C on , 
and F on , see Fig. 12.31 Nevertheless, we have found offsite cases where these returns do occur in the 
second finite gap, e.g. at a « 0.861 and T+ = 1. 

Homotopy between Onsite and Offsite GSs.. As one can see in Fig. l2.3l and T2.4i SGS families often 
intersect the vertical line T_ = r+, which describes a no- interface medium. The numerical results 
show that at each of these intersections the corresponding solution is a symmetric GS. Clearly, for 
(|2.4[) with xo = (used in SGS bifurcations from onsite GSs) if a GS is centered at x = nd, n € Z, it 
is an onsite GS and if it is centered at x = (2n + 1)|, u6Z. it is an offsite GS. In the case xq = i 
(used in SGS bifurcations from offsite GSs) if the center of a GS is at x = nd or x = (2n + it 
is an offsite or onsite GS respectively. Fig. 12.81 plots the GS solutions located on the vertical line 
r_ = r + for the families C on and F Q ff. It also shows the center of mass for all the SGSs along the 
families clearly visualizing that the SGS is no always shifted to the more focusing medium. 

Although in Fig. 12.81 the function xcm{4>) appears to be monotonic in the arclength r, this is not 
always the case; not even for open SGS families. Fig. 12.91 shows the family E on from Fig. 12.31 (e) for 



Surface Gap Solitons; Stability via the Numerical Evans Function 



<icr 4 (a) r -1 




(b) r.1 







G ff 1 




off 




2 r -1 



Fig. 2.4. SGS families bifurcating from the off site GSs A g— H g in Fig. \2.2\ Note that the curves were computed 
up to |r_| = 30 or \Efj,(rp)\ = 100, whichever occurred first, and the behavior is qualitatively the same as in the displayed 
windows. 
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Fig. 2.5. Profiles of SGSs in the C g family, (a) the SGS family in the (r_ , (</>)) plane with 9 points labeled; 
(b) and (c) profiles of the SGSs labeled in (a). 



which four intersections of T_ = T + occur but three of these are GSs centered at x = d/2 = 5. Note 
also that these three GSs are all different. 

The numerics suggest that each shift of the SGS center from one extremum of V to another 
is accompanied by a fold. Moreover, for spatially relatively broad SGSs. i.e., when \x lies near the 
corresponding bifurcation edge, all folds seem to correspond to such shifts, see Fig. I2.8l and [2.10I The 
same occurs also for the cases A on , F on , A a g and C ff. 
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r xx 



Fig. 2.6. Translated bifurcation GSs within SGS families. Left: Families of SGSs in the case C on and F on at the 
top and bottom resp. Middle: GS at the (first) return of the full line in the left plot to (T+, En(<j>Gs))- These GSs, 
denoted GSi are centered at x = d. Right: GS at the (first) return of the dashed line in the left plot to (T+, E^(ij>Qs)) ■ 
These GSs, denoted GS—±, are centered at x = — d. 




Fig. 2.7. Two GSs centered at x = within the E g family, (a) (r_, i? M (</>)) curve of the SGS family; (b) center 
of mass; (c) profiles of the two solutions on the Y— = T+ line. Both are GSs centered at x = 0. 



Asymptotics near Bifurcation Gap Edges.. It was shown in [13] that in the vicinity of spectral 
gap edges (J, = s n the interval of existence T_ e (r + — a(/i), T + + b(fi)) with a(fj), b(n) > for SGSs 
localized near x = collapses as fi — > s n . In detail a(/x), b(fi) — > 0+ as /i — > s n — for T + > and as 
(j, — > s„+ for r + < 0. In practice, since folds typically correspond to transitions of SGS centers to 
the next extremum location of V, this means that the horizontal distance |T_ — T + \ of any GS point 
along the SGS family to the nearest fold converges to as the edge is approached. Here we show that 
the same holds for SGSs localized around any point x c 6 M. Analogously to Sec. 3.3 in [13] we use 
the asymptotic representation 



(j, = s n + e 2 n + 0(e 4 ) 
4>(x) = eA(X)q n (x) + e 2 A'(X)q n (x) + e z ^{x,X) + 0(e 4 ) 



(2.5) 
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Fig. 2.8. On- and off site GSs within an SGS family. Top and bottom: SGS families for the case C on and F g 
resp. Left: total energy (</>); middle: center of mass xcm(<P)> right: profiles of GSs marked on the left and in the 
middle. For the case Con (where xq = 0) profiles GS n with integer indices n £ Z are onsite GSs while those with 
indices n/2,n G Z are offsite GSs. For F a g (where xo = 5) this is vice versa. The dashed vertical lines correspond to 
a no interface medium (T_ = ). 
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Fig. 2.9. SGS family through E on ; description as in Fig. \2.8\ Notice the non-monotonous dependence of the 
center of mass on the arclength in contrast to the cases in Fig. \2.8\ 



with X = e(x — x c ), < e « 1. Substituting this two-scale expansion in (|1.5[) . we get at O(e) and 
0(e 2 ) 



(L - s n )q n = and (L - s n )q n = 2q' n 
respectively, which both have periodic solutions. At 0(e 3 ) we have 

(L - s„)0 (3) = nA n q n + A"(q n + 2q' n ) + T(X + ex c )A 3 q 3 n , 



(2.6) 
(2.7) 



where we have used the fact that due to the form T(x) = T + X[o,oo)( x ) + T- Xf-oo.ojl 2 ') it is T(x) 
T(ex) = T(X + ex c ). 
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Fig. 2.10. Behavior near gap edges. SGS family bifurcation from the onsite GS at fi 0.28242 (near the minimum 
of the first band) with = 1. Left: total energy E^^tf)); middle: zoomed-in parts of the plot on the left; right: center 
of mass. Notice the large number of folds near the GS bifurcation point (r_ , (</>)) = E^^gs)) ■ 



To ensure the existence of a periodic solution cjP^ , the Fredholm alternative needs to be applied 
to (|2.7[) resulting in 

fL4 + vA" + pT{X + ex c )A 3 = (2.8) 

with v and fi defined in (|2.3|) . 

Due to Lemma fTTTl cquation (|2.8p has non-trivial localized solutions only for T + = T_. For the SGS 
families this means that as \x approaches the corresponding spectral edge s n , they lie inside narrower 
and narrower vertical slabs T_ e (r + — a(/i), T + + 6(/i)) with a(fi), — > as fi approaches s n . In 
other words the horizontal distance |T_ — T + | of every fold to the GS bifurcation point converges to 
as (i approaches s n . The convergence rate is, of course, expected to decrease with \n\, where n is the 
signed ordinal of the fold with respect of the bifurcation GS, since the effect of the interface decreases 
with increasing distance between the solution center and the interface location x — 0. In fact, one 
can conjecture an exponentially fast decrease of this rate as \n\ increases due to the exponential decay 
rate of SGS tails. In other words, we expect that for a fixed // with — s n \ << 1 the distance of the 
n-th fold grows exponentially with |n|, cf. the right panel of Figure l2".10l 

For GSs centered at x — the first fold in each T_ direction was followed numerically in fi in [13| . 
Although we have not followed other folds systematically, our numerical results are in agreement with 
this asymptotic statement. Fig. 12.101 presents a part of an SGS family bifurcating from the onsite 
GS at ii w 0.28242, i.e., close to the edge si 0.28317. The first 4 folds in both directions |n| < 4 
occur for |F_ — T + | < 0.13. The part of the family between the folds n = —4 and n = 4 consists of 
solutions with centers of mass in (—20, 20). It thus contains also solutions centered relatively far from 
the origin. For /i w 0.28278, i.e., even closer to si, solutions between the first 18 folds (|n| < 9) have 
|r_ — r+| < 0.1 and centers of mass in (—50,50). The plot of the curve (T - , E ^((f))) is even more 
complicated and less readable than in Fig. 12.101 

3. Linear Stability of SGSs. One of the most important properties of the above computed 
SGSs of (jl.ip is their stability with respect to perturbations of initial data. Only stable SGSs can be 
viewed as physically relevant states of the system. Orbital stability of bound states e -1AIs *</> s (a;) of 
with <p s positive and [i s in the semi-infinite gap can be checked using the Vakhitov-Kolokolov criterion 
on the sign of ^^^^ , where P(<f>(-,[i s )) = ||</>(-,A i s)|||2( R ), and on the kernel and the negative 
part of the spectrum of the operator L + (defined below in the proof of Lemma 13. ip [HJ [HJ HSj • In 
detail, orbital stability is satisfied if dp ^(-.M 3 )) < q an( j ^ + ^as no zero eigenvalues and the number 
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of its negative eigenvalues is 1. These conditions have to be in general checked numerically and, in 
particular, the spectral condition on L + is not completely trivial to check. Moreover, this test applies 
only in the semi-infinite gap. In order to be able to make stability statements even about SGSs in 
finite gaps, we choose to study linear stability by inspecting directly the spectrum of the linearized 
operator via the numerical Evans function method. 

GivenaSGS e~ 1,J,st (f> s (x), we consider the evolution of a perturbed solution ip (x, t) = e~ lflst ((f> s (x)+ 
q(x,t)) with q(x,0) small. Linearizing (|1.1[) in q and using that (fx s ,4>s) satisfies (|1.5[) . we obtain 

ld t q + H s q + d 2 x q-V(x)q + T{x)4> 2 a (x){q + 2q) = 0. (3.1) 

As a system for (J7i, U 2 ) ■= (q, q) this becomes 

Ui\ . fUi 



u 2 j= h {u 2 ) 

with 

- /i., + V(x) - 2T(x)tf -n*)4> 2 s \ (o 

y r(x)tf dl + ^-v(x) + 2T(x)^) ^- 6 > 

As ()3.2[) is t— autonomous, the problem is separable and we can let (Ui(x, t),U2(x, t)) = e xt (u\{x), 112 (x)) 
so that we obtain the spectral problem 

*(:)=>-&)• <3 - 4 > 

Clearly, linear stability of e~ ltlst <f) s {x) is determined by the real part of <t(L). 

3.1. The Spectrum of L. Although cr(L) consists of the essential spectrum a ess (L) and of 
eigenvalues CTdi SC (L), only eigenvalues dictate stability properties in this case because a ess (h) C DR. 
This follows from Weyl's theorem, see Sec. XIII. 4 of [42], which ensures that due to the exponential 
decay of <j) s , the terms containing <p s are a relatively compact perturbation of L and the essential 
spectrum of L equals that of i ^ L ^ ^j- As this operator is skew-adjoint, its spectrum is 

imaginary. Moreover, we also have that 

er css (L) = i(a(L) - fi s ) U —i(a(L) - fj, s ), 

where the first bands of cr(L) are detailed at the beginning of Section |2~T1 Fig. 13.11 plots <7 eS s(L) for 
the potential (|2.4[) schematically. 

Although CTdi S c(L) will be determined numerically, the following Lemma gives an upper bound on 
Re((7disc(L)). For the dynamics of (|3.2[) this is an upper bound on the growth rate of the perturbation. 
Note, however, that the numerical method for studying crdisc(L) will not make use of this result. 

Lemma 3.1. Eigenvalues A £ a disc {l^) satisfy |iZe(A)| < | [r| 1 00 1 1 0s 1 1 So - 

Proof. With the new set of variables a :— u± + u 2 and b := —i(ui — u 2 ) equation (|3.4j) becomes 

\( a \-( ° M ( a \ where L = ~ d * + V{x) ~ ^ ~ (3 5 ) 
X {b)-{-L + 0){b) whcrC L + =-dl + V{x)-^-ZY(x)4>l{x). (3 - 5) 

After multiplication of the first equation in (|3.5[) by a and of the complex conjugate of the second 
equation by 6, integration over R and addition we obtain 

A||a||| 2 + A||6||| 2 = / -b"{x)a{x) + a"{x)b{x) + 2T{x)4> 2 s {x)a{x)b{x)dx. 
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Fig. 3.1. The essential spectrum o/L for V(x) = sin 2 (7ra;/10) . 



Integration by parts leads to 



|A||a 



| 2 , 2 +A| 



- I 

\l 2 \ 



T(x)<f>l (x)a(x)b(x)dx 



< 2iin 



I 2 1 

I GO I 



a \\m 



\L 2 , 



where the inequality follows by Cauchy-Schwarz. Finally, for Re(A) we obtain |Re(A)| < 1 1 T 1 1 oo 1 1 <^>s 1 1 
using the inequality 2a(3 < a 2 + f3 2 for a, (3 G R. □ 

Earlier works apply a discretization method to equation (|3.1[) and compute the eigenvalues of 
the corresponding system directly, see e.g. |47j . where GSs were studied. Yet, as pointed out in 
[4"T1 [T2"] . this may be highly computationally costly and typically leads to spurious eigenvalues of 
(|3.4[) . Identifying those spurious values is a nontrivial and absolutely necessary task when studying 
stability of solitary waves. The Evans function method was shown in |12| I41j to avoid this problem. 



3.2. The Evans Function Method. The Evans function £(X),£ : C — > C is a generalization 
of the characteristic polynomial for a linear operator. Its zeros coincide with isolated eigenvalues of 
the operator. The multiplicity of the zero equals the algebraic multiplicity of the eigenvalue. The 
Evans function is defined and analytic away from the essential spectrum of L. It was introduced first 
by J. Evans in his study of stability of nerve impulses [17] . For analysis of the Evans function see, 
e.g. 0110]. It has been used for both analytical H3 HQ [25] and numerical [H E] E] SU [19] studies 
of stability of traveling and standing solitary waves. For problems with exponential dichotomy the 
Evans function is a determinant of a matrix generated by bases of the stable and unstable manifolds 
corresponding to the trivial solution u = of (|3.4p . The stable and unstable manifolds consist of 
solutions that decay exponentially fast as x — > oo and x — > — oo, respectively. When this determinant 
vanishes at a given A, the manifolds are linearly dependent, which means that a solution exists which 
decays exponentially for both x — > oo and x — > — oo, implying A € (Tdisc(L). 

In |41j the numerical Evans function method was used in combination with the variation of 
constants to study linear stability of GSs of (|1.1[) (with T + = The authors, however, did not 

search systematically for eigenvalues over the whole right half plane. Instead they concentrate on real 
eigenvalues and on eigenvalues bifurcating from the edges of er 0SS (L). We show, in fact, that for a 
numerically stable and accurate evaluation of £ (A) for |A| large a change of variables and the use of 
exterior algebra are necessary. With the help of these tools, which were not applied in [41], we carry 
out winding number computations of £ (A) with a contour 7 = 5 + iR, < 5 << 1, parallel to the 
imaginary axis in order to determine the number of eigenvalues in the whole half plane to the right 
of 7. Note that thanks to an asymptotic behavior of £ (A) for |A| — > 00, we will be able to use a finite 
contour. 
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3.2.1. Construction of the Evans Function for (|3 .4[) . As advertised above, the Evans func- 
tion is constructed using the stable and unstable manifolds of the trivial solution u = of (|3.4p . For 
|z| — ^ oo the terms in (|3.4[) that are proportional to <p s {x) vanish due to the exponential decay of <f> s 
and one obtains the uncoupled Hill's equations 



u%.u\ + V(x)u\ = n+ui 



(3.6a) 



— d^.u 2 + V(x)u 2 = H-U2, (3.6b) 
where fj,± := /i s ± iA. If yU+ ^ ^{L), the fundamental system of (|3.6ap consists of the Bloch waves 

< + 2 (z) = Pt%(x)e ±ik + X , with pfi(x + d)= p$(x), Im(fc+) > 0, (3.7) 
which decay exponentially fast asi-> ±oo respectively Similarly for ^ tr(L) the Bloch waves 

^1,2 (&) = Pi72( x ) e±ik ~ X with Pia( x + d ) = Pi.ifr), Im(fc-) > 0, (3.8) 

that build the fundamental system of (|3.6b[) . decay exponentially fast as x — > ±oo respectively. We 
choose the following normalization of the Bloch waves in (|3.7[) . (|3.8[) 



< + 2 (0) = l, <2(0) = 1- 



(3.9) 



Note that for even potentials V this normalization is possible since n± do not lie in a(L) and p\^ 
arc thus not Dirichlct cigcnfunctions of () 1 . 8[) (with k = k+ and k = fc_ for p^ + 2 and p^ ~ 2 respectively) 
so that Pi 2 (0) 7^ 0. For an even V a Dirichlet eigenfunction p(x) would allow the extension of the 
solution ip(x) = p(x)e lkx via symmetry from the half line on which -0 is bounded to K, producing so 
a bounded solution, which is impossible since fj,± ^ o~(L). 

For x — > oo the stable manifold of (|3.4p is thus generated by 







and for x 



-oo the unstable manifold of (|3.4p is generated by 







(3.10) 



(3.11) 



The four vector functions in (|3.10p and (|3.1ip are asymptotic reductions of the fundamental solution 
set of p.4p though each of the pairs is for a different asymptotic region. Note that we use the 
normalization (|3.9p also as a normalization of the fundamental solutions of (13.41) . 
Let us rewrite (|3.4p as the first order system 



Av, 



A 





V{x)~2T{x)4> 2 s (x)-^ + 




V 



-T(x)tf(x) 



1 

-T(x)tf{x) 



V{x)-2T{x)4> 2 s {x)- ^ 




1 



(3.12) 



for i> := (ux, u^, u%, u' 2 ). The fundamental system of (|3.12p is denoted by {v x , v 2 , , v%}- The labels 
are assigned according to the above asymptotic behavior as follows 



v 1 (x) 



( \ 






v 2 (x) 



as x 



and 



(3.13a) 
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Note that A is a parameter of v 12 . The Evans function £ is then defined as 
£(X) = dct(vi (a;*; A), v^ix*; A), (x*; A), v^ix*; A)). 



(3.13b) 



(3.14) 



As £ is constant in x*, its choice is arbitrary and for our purposes of stability of SGSs we choose 
= 0, i.e., at the interface location, which avoids the need of numerical integrations of (|3.12[) across 
the interface. 

The basic (naive) idea of how to evaluate £ (A) is to solve Q3.12p for v^ 2 up to x = x* with the 
initial conditions given by (|3.13ap evaluated at x = —L^ << — 1 and for v^ 2 up to x = x* with the 
initial conditions given by (|3.13b[) evaluated at x — L^. In Section we discuss difficulties of this 
approach. 

A-priori information about the behavior of £ includes the upper bound on the real part of its 
zeros given by Lemma l3.H Besides, we have the following symmetry of £ and asymptotic behavior of 
£(X) for |A| oo. 

Lemma 3.2. The Evans function £ in (|3.14p satisfies £(X) = £(X). 

Proof. Complex conjugation of (|3.12[) reveals that v(x; A) and Pv(x; A) satisfy the same differential 



equation u' = A(X)u, where P is the permutation matrix P 



10 
1 
10 
10 



Because ipi + < " X \x) 



ipi (x) and ip2 + (x) — ^2 the asymptotic data (|3.13|) arc related via complex conjugation 

and the permutation P in the following way: v^(x, A) ~ Pi^(a;,A) as x — > — oo and Vj~(x,X) ~ 
Pv 2 {x, A) as x — > oo. We thus have 



£ ( A) = dot (v-f (.t*; A), (x*; A), (x* ; A) , v£(x*; A)) 

= dct (P(ti^(x*; A), uf (x*; A), v 2 (x*; A), Vi(x*; A))) 

= det (P(i>f (a;*; X),v 2 (x^; X),vf (x*;X),v£ (x*;X))P 2 ) = £{X) 



where Pi := 



10 
10 
1 
10 



and where the product rule for determinants together with the fact that 

det(P) = dct(P 2 ) = 1 were used. □ 

Lemma 3.3. The renormalized Evans function 



E(X) := 



-1 



4A + 1 



£{\) 



(3.15) 



is analytic for Re(X) > and satisfies E(X) — > 1 as |A| — > oo for arg(X) G (— tt/2, 7r/2). 
Proof. The analyticity follows from analyticity of 4 ~_ | l ; 1 and £ (A). 

For |A| large all terms except those proportional to A become negligible in the second and fourth 
equations in (|3.12j) leading to 



Av with A 



10 
-iA 
1 
iA 



(3.16) 



The 4 eigenvalues of A are ±\/—iX and ±viA with eigenvectors (1, ±\A-iA, 0, 0) and (0, 0, 1, ±viA) 
respectively. The fundamental solution set of (|3.16p can thus be chosen as 



-iAa~ 



1 






, e 



/iXx 






1 



, e 



-V-iAa: / -^-iA 





, e 



— V i\x 






1 

-\/iA 



(3.17) 
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With the convention that the complex argument of the '+' square root function y/z lies in (— it/2, tt/2], 
i.e., arg( v / i) 6 (— tt/2, tt/2], we get by comparing with the asymptotic behavior in (|3.13a[) . (|3. 13b|) 
that the ordering of the vectors (|3.17[) agrees with the one in £(X). Also the normalization agrees 
with the one implied by (|3.9[) . Hence the vectors in (|3.17[) can be used directly to evaluate £(X) for 
|A| large. A calculation of the determinant of the vectors in (|3.17[) for |A| — > oo yields £(X) ~ — 4A if 
arg(A) € (—tt/2, tt/2) and £(X) ~ 4A if arg(A) € (tt/2, 37r/2). The imaginary axis is excluded to avoid 
the essential spectrum where our Evans function has not been defined. The renormalized function 
E(X) thus converges to 1 for arg(A) G (—tt/2, tt/2). □ 

3.3. The Winding Number Method and a Choice of the Contour. The argument prin- 
ciple |38j is a standard tool for counting the number of zeros of a complex function within a region of 
its analyticity. Denoting by n(7, z) the winding number [38j of a closed curve 7 C C with respect to 
the point z € C \ 7, the argument principle states that if E is an analytic function inside 7 and has 
no zeros on 7, then 



where N is the number of zeros of E, counting multiplicities, inside 7. 

Using (|3.18|) we can, in principle, determine the number of eigenvalues of L in the right half 
complex plane by using a contour 7 that encloses it. £ is, however, not analytic on er css (L) C 1R so 
that 7 has to stay a positive distance away from the imaginary axis. Although it should be possible 
to extend the Evans function analytically across the essential spectrum by generalizing the analysis 
of [331 HH for our periodic coefficient case, the use of the imaginary axis as a contour for the 
argument principle would still be unpractical due to the expected occurrence of zeros of the extended 
Evans function within the essential spectrum. This expectation is based on the analysis of constant 
coefficient problems, for instance, in [24], where zeros in the essential spectrum arc found and their 
bifurcation from spectral edges studied. 

We take 7 j| iffi with 5 := dist(7, UR) > 0. In most of the numerical computations we use S = 0.005. 
Due to the <5-gap between 7 and the imaginary axis, we are thus able to detect only eigenvalues with 
real part greater than <5, i.e., only 'substantial' instabilities. Using the symmetry from Lemma 13.21 and 
the asymptotic result from Lemma 13.31 the curve 7 can in practice be replaced by 



where H > is chosen for a given <p s so that near Im(7) = H the Evans function E(X) is clearly 
converging to 1. In most of our computations H = 20 is sufficient. In order to numerically evaluate 
the winding number, we have written a simple Matlab script, see Appendix [A"l 

It is possible that E has a zero on 7, which would prevent the applicability of the argument 
principle. In numerics this will, however, generically result in a small value of E(X) so that the effect 
is the same as a zero near 7. If E attains numerically the value somewhere along 7, the contour can 
be locally deformed. This case does, however, not occur in our simulations. 

Reliable computations of n(E(j),0) require a sufficiently fine discretization of 7. Needlessly fine 
discretization, on the other hand, lead to intolerably long computation times. A possible solution is 
the adaptive discretization suggested in Sec. 3.3 of [4] with the main idea being to ensure that the 
difference in the complex argument &rg(E(X)) for neighboring points Ai 2 G 7 is near a desired value. 

3.4. Numerically Stable Evaluation of E(X). Attempting a straightforward evaluation of 
E(X) as defined via Q3.12p - p,15|) . one encounters several issues with numerical stability and accuracy. 
Firstly, the evaluation of the initial conditions (|3.13[) at some x = ±Loo, >> 1 requires computa- 
tion of exponentially growing/decaying Bloch functions, where the growth/decay rate is very strong 
for certain parameter values. Secondly, the solutions of (|3.12[) also feature exponential growth/decay. 
In both of these cases parameter dependent transformations can be used to remove the exponential 



n(E( 1 ),0) = N, 



(3.18) 



7 = <5 + i[0,iJ], 



(3.19) 
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growth/decay. Next, more importantly, the system (|3.12l) suffers from stiffness when the growth/de- 
cay rates of the two vectors in the manifold of interest (stable or unstable) are largely different. This 
can be overcome by a reformulation of (|3.12[) in the exterior algebra. As explained below, the use of 
Grassmanian preserving ODE integrators is then vital. In the following we describe these problems 
and their solutions in detail. 

3.4.1. Computation of the Bloch Functions in (|3TT3|) . Due to the form ([3?f|>. (|3~8l) the Bloch 
functions need to be, of course, computed only on the period x £ [0,d]. In (|3.13p their values at 
x = iLoo, Loo >> lj which are used as initial data for (|3.12p . are then obtained straightforwardly due 
to periodicity of p 1 ^ \ . Two commonly used methods are available to compute the Bloch functions. The 
first possibility is to solve (|1.8p with u>(k) = /Lt± fixed as a quadratic eigenvalue problem in (k,p ,1± ) 
on x G [0, d] with periodic boundary conditions. To avoid the need for solving quadratic eigenvalue 
problems, we opt for the second traditional method, which is Floquet theory [33JQ3], in which (|3.6a[) 
and (|3.6b[) arc solved on [0, d] as initial value problems and k± are determined from eigenvalues of the 
monodromy matrix, i.e., from Floquet multipliers. Even on the bounded domain [0, d] the exponential 
growth of i/i[2, however, leads to a strong enhancement of error when |Im(A;±)| is large. As shown 
below, this occurs, in particular, for \/jl±\ large when arg(/x±) is not close to tt or — tt. 

For \n\ large the equation —d^w + V(x)w = fj,w reduces effectively to 

—d 2 .w = fiw 

with the fundamental system {e ± ^ = ^ x }. With the definition of the '+' square root y/z, z g C such 
that &Tg(y/z) G (— 7r/2, it/2], the solution e^~~^ x grows fast when is large and arg(/x) is not close 
to ±7r. In (|3.6|) and with the choice of the contour along the upper half of the imaginary axis as in 
(|3.19j) this happens for /i_ (= fi s — iA) when |A| is large. In (|3.8[) we thus get a large Im(fc_). For 
Li + (= fi s + iA) the contour ()3.19|) leads to arg(/i + ) near tt for \(j, + \ large so that no strong growth 
occurs in i- e -j Im(fc_) is close to 0. 

To achieve sufficient accuracy of ip^ ~ 2 , we use a change of variables that removes the exponential 
growth. Rewriting (|3.6b[) as a first order system for v := (u2,u' 2 ) T , we have 

v' = Bv, B:=( v{x) ^^ j). (3.20) 

The Bloch functions of (|3.20[) arc 

v^(x) = «-(*), «-)'(*)) T , v&(x) = {^-{x)M-)'{x)) T , 
which due to (|3.8[) have the form 

V W( X ) = qW(x)e ik - x , v {2 \x) = q {2) \x)e-' lk - x , (3.21) 

where 

qU(x) = ((t-vr^JS f— ( ^\ and q™(x) = ( , „_ w f ~ 0=) \ ( 3 . 2 2) 

To remove the exponential growth, we introduce v^\x) := {x)e~^~^~ x such that 



v' = {B - y/^fjZl)v. (3.23) 
Via standard Floquet theory [14] the Bloch functions 



(3.24) 
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of (|3.23[) can be computed by solving (|3.23[) on x 6 [0, d] and finding the eigenvalues of the monodromy 
matrix, i.e. the characteristic multipliers p\ i2 = e lkl < 2d . As a 'side effect' the change of variables v — > v 
has doubled the decay rate in which possibly results in extremely small values of near 

x = d such that for large it is numerically impossible to recover accurate information about 
We, therefore, use (|3.23l) to compute only v^ 2 \ The corresponding is identified as the solution 
belonging to the larger characteristic multiplier, i.e., \p^\ > \p^\, and k 2 is then obtained from 
k 2 = — ilog(p^)/d. The periodic part is computed from (|3.24j) . Finally, k_ is given by 

fc_ = iy/—fi- - k 2 . 

The missing Bloch function or cquivalently ip^~~ , is readily available via symmetry. For 
V'li = Pi 2 (x)e ±ik ~ x equation ()3.6bj) becomes 

- (d x + ifc_ ) 2 P t- + V(x)pt- = fi-pf- , (3.25a) 

- (d x - ik-fpz 1 - + V(x)p^- = pi-pl~ (3.25b) 

with periodic boundary conditions on a; £ [0,d]. The substitution fc_ — > — fc_ and x — > —x in the 
second equation reveals that 

Note that p^~ is unique up to a normalization, i.e., p_ is a simple eigenvalue of (|3.25a| since otherwise 
for fc_ =/: one could generate via the above symmetry at least four linearly independent solutions 
of (|3.6b[) . The case fc_ = yields a periodic Bloch function and can thus only occur for /i_ E &{L) 
[33l[T4], which is not our case since Im(/x_) = — Re(A) ^ 

3.4.2. Dealing with the Stiffness of (|3.12p . In Sec. (|3.2.1|) the Evans function is defined 
using the bases {wf , v 2 } and {vf, v£} of the unstable and stable manifolds respectively of the trivial 
solution of p,12[) . In practice i>-j~ 2 need to be computed from x = —Loo << —1 up to x = x* and v^~ 2 
from x = Loo >> 1 down to x — x„. In the unstable manifold and v 2 have different growth rates. 
For large where the terms proportional to <j>1{x) are negligible, the rates are e lva ( k +) x and e Iui ( k -) x 
respectively. When lm(k + ) and Im(fc_) differ strongly, the problem is stiff and the fast growing 
solution dominates the dynamics and in a numerical simulation the error is enhanced in the direction 
of this more unstable mode so that the solution growing more slowly cannot be accurately computed. 
The functions and v 2 thus do not remain linearly independent in the numerical x— evolution 
although they are linearly independent initially, see p.!3a[) . An analogous situation occurs for 
and v 2 of the stable manifold. We saw in Sec. 13.4.11 that for the selected vertical contour this large 
difference of growth rates occurs, in particular, for Im(A) large, where Im(fc_) is large and Im(/c + ) is 
close to zero. 

A traditional approach to a problem with linear dependence is to apply orthogonalization. Or- 
thogonality for all x is, however, not generally satisfied by solutions of (|3.12j) . (|3.13j) . and as explained 
in [3"1 1 1 2 j . orthogonalization may destroy the analyticity of E = E{\). This is inadmissible when we 
intend to employ the argument principle. An alternative is the use of exterior algebra as first proposed 
in [37] for solving stiff linear ODE systems. More recently it has been used in spectral problems, for 
example, in [4j [3j [12l US] . For theoretical background on exterior algebra see, e.g., [20l IFJl [50] . 

Let {ei, . . . , e„} be an orthonormal basis for C™. We first define the wedge (exterior) product via 
the properties of distributivity and associativity 

e l A (aej + e k ) =aei A e 3 + e.; A e k 
(aei + ej) A e k —aei A ej + ej A e k 
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and antisymmetry 

e, A ej — —ej A ej. 

Then the distinct and nonzero members of {e,j A • • • A e.; fc : 1 < i\ < ■ ■ ■ < ik < n, 1 < k < n} form a 
basis of the vector space of all fc-forms A ... A cS- k > with at 3 ' <e C n ,j = 1, . . . , k, which is generally 
denoted by /\ k (C n ) and has dimension (^). Let = J27=i ^ij e i with j = 1, • • • ,n, and B := 
Then (see [SO], Proposition 3.3) 

A ... A 6 (n) = det B e 1 A . . . A e„. (3.26) 
For example, for a, b <G C™, a = Y^7=i a i e ii b = S™=i one ^ nus nas 

n— 1 n 

a Ab = (a^j — djbi)ei A e^-. 

Z— 1 j — 2+1 

Furthermore, for m S /\ J (C") and u e /\ fe (C") one obtains mAd 6 /\ 3+fe (C"). In particular, for 
j = n — k we get iiAug /\™(C"). In this case (see [3], Section 5) 

«Au = (S, Eu)ei A ... A e„, (3.27) 

where (■, ■) is the standard inner product on eft) with conjugation of the first argument and S is an 
orthogonal matrix. 

Applying (|3.26p to our Evans function problem with k = 2, n = 4, we get 

E(\)e\ h ... hen = — -ur(x„; A) A v^ix*; A) A v^(x*; A) A (a;*; A) 

4A -)- 1 

= ^ I V-(x.;A)AV+(x,;A) ) 

where := o^A^. Because £ /\ 2 (C 4 ), the rule p.27|) applies and 

E(X) = ^- Y (W,EV + ). (3.28) 

Indeed, -E(A) = if and only if v± ,vf, and v£ are linearly dependent (see e.g. [SO], Lemma 3.4). 
In [3] it is shown that the form of £ for this case is 

( 1 \ 
-1 1 
o'o'ioo'o 1 (3-29) 
-1 / 
1 0/ 

provided the basis for /\ 2 (C 4 ) has been selected as e\ A e 2 , e\ A e^, e\ A e 4 , e 2 A e 3 , e 2 A e 4 , e3 A e 4 . 
Evolution equations for are obtained by linearity of the wedge product 

V^' = vf A + vf A vf = Avf Avf +vf A Avf. 

As described in Section 2 of [3], denoting by ajk the (j, fc) entry of A, this induces the matrix 

' an+a 2 2 a-23 a 24 — a-13 -Ol4 

0.32 011+033 «34 0-12 -«14 

4(2) _ I «42 a 43 an+a 44 a i2 a i3 

— 031 "221 022+033 ^34 — «24 

— a 4 i a 2 i 043 022+044 Q23 
— 041 a 3 i — a 42 a 32 033+044 
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such that 

y±' = A (2)y±, (3.30) 

An alternative to solving p. 121) for each of v± and from x — — to x — x* and and v^ - 
from x = Lqo to x = 1, is thus to solve (|3.30p for V~ from x = —Loo to x = x* and for V + from 
x = Loo to x = x*. The 2-dimensional manifolds of (|3.12p are represented by lines in p.30[) and in 
this way we avoid the linear dependence problem discussed above. The initial conditions for (|3.30p at 
x = ±Loo are determined by the wedge product of respective vectors in (|3. 131) . 

Decomposability and preservation of the Grassmanian. . It is not true that every k— form can be 
written as a single wedge product of k elements of C™. In other words not every k— form represents a 
fc— dimensional subspace of C™. If this is satisfied, i.e., if U G /\ fe (C n ) can be written as 

U = A ... A 

with u^l G C n ,j = l,...,k, then U is decomposable [2171 IrJl 13] . 

In our setting are decomposable initially at x = ±Loo by construction, but in a numerical 
integration of (|3.30p this property may not be preserved so that at x — the computed solution no 
longer represents the 2-dimcnsional stable or unstable manifold respectively. Restricting to n = 4, k = 
2, the following result holds. A 2-form U G /\ 2 (C 4 ) is decomposable if and only if U A U = [20] . For 
the coefficients it means 

!{U) := hu, EU) = UiUe - U 2 U 5 + U 3 U 4 = 0. 

The quadratic quantity I{U) is called the Grassmanian. It is a strong invariant of (|3.30p if T,A^ + 
(T,A^) T = 0, see [3], which holds in our case, where A^> is constructed from A defined in (|3.12l) . 

A number of methods are available that preserve strong quadratic invariants |22j . As discussed 
in Section l3.4.3[ we, however, solve (|3.30p after the change of variables 

y± ._ e ±™ v ± ( 3 31 ) 

with v = Im(fc+)+Im(fc_). With (EOH we have I(V ± ) = 7(y ± )e ±2ra and ^(/(V^)) = £(I(V ± ))± 
2vI(V ± ) = ±2vI(V ± ) so that / is only a weak invariant of the equation for V^, namely equation 
(|3.32p . This means that if I(F ± (a;o)) = at some initial point x , then 7(V r± (x)) = Va; G K. Weak 
invariants are generally impossible to numerically preserve exactly. Approximate preservation can be 
achieved by projection methods |22| . In our case the numerical solution is projected onto the manifold 
M := {U G C 6 : I(U) = 0}. We solve (|3T32|) via the standard 4-5th order Rungc-Kutta method of 
Matlab (ode45) and execute the projection at regular x— intervals, usually of length Sp = 0.2 (i.e. the 
frequency of the projections within one period of the potential V is d/Sp = 50), see Sec. 13.5.11 for 
more details on the choice of Sp. 

Renaming the solution vector of (|3.32p to w, the projection constitutes the following optimization 
problem. Given a numerical approximation w (at some point x) , find w € M such that 

— w>|| = min \\U — w\\, 
ueM 

where || ■ || denotes the Euclidean norm in C 6 . Defining the Lagrange function \\w — w|| 2 /2 — £/(«;), 
we get the following nonlinear problem for the Lagrange multiplier £ 

I(w + (I'(w) T ) = 0, 

which we solve via Newton's iteration 



£(n+l) _ ^(n) 
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where 

■^(I(w + (I'(ti) T ))(() = w\ + 2Cu>iu> 6 + w\ + w\ - 2Cw 2 w 5 + w\ + w\ + 2Cw 3 ^4 + w\. 

In Section 13.5.11 we show that the accuracy of the Evans function increases with increasing pro- 
jection frequency d/Sp. 

3.4.3. Removing the Exponential Growth in (|3.30[) . Based on the asymptotic behav- 
ior p,13[) of the solutions vf 2 an d linearity of the wedge product, we get that V~ behaves like 
e (imO + )+im(fc_)):r j? or x \ aT g C negative and V + behaves like e -( Ini ( k +)+ Iml , k -)) x for x large positive. 
This exponential growth in the backward and forward x direction for respectively leads to diffi- 
culties with numerical accuracy when Im(fc + ) + Im(/c_) is large. Clearly, for values of A close to zeros 
of E the growth is eliminated in the respective opposite asymptotic regions (at x » 1 for V~ and at 
x « —1 for V + ) as the stable and unstable manifolds are close to being linearly dependent. In the 
intermediate part of the domain (e.g. near x = x*) the growth is possibly largely limited. At other 
values of A the growth may, however, remain large for all x. We, therefore, perform the numerics of 
(|3.30[) under the transformation 

y± ._ e ±(Im(fe+)+Im(fc_))xy±_ 

The new variables satisfy 

V*' = (a^ ± (Im(fc+) + Im(fc_))/) V ± . (3.32) 

Linearity of the wedge product yields 

E(X) = ^TT V ~( X *^) A V+(x*;\) = ^L_V-(x,; A) A V+(x„; A). 

3.5. Numerical Results on Stability of SGSs. This section presents results of Evans function 
computations for the spectral problem corresponding to linear stability of the SGS families in Figs. 
12.31 and 12.41 In several examples the NLS (|1.1[) is then also evolved numerically in time with the 
selected SGS profile as the initial condition in order to demonstrate the stable or unstable behavior. 

3.5.1. Accuracy of the Numerical Evans Function. The problem at hand offers a natural 
test, namely E(0) = 0, which holds because is an eigenvalue of L and ^ cr e ss(L). This kernel is 
caused by the phase invariance of (jl.ip . In the numerics we use |-E(0)| as a measure of accuracy of 
the numerical Evans function. There is a number of factors that may affect accuracy: (i) error in 
computing the Bloch functions of p.23[) ; (ii) numerical error in solving p.30[) for the two manifolds; 
(iii) size 2Loo of the computational domain for (|3.30[) ; (iv) norm of the residual in computing the 
SGS profile S ; (v) h and p in the FEM computation of the SGS profile; and (vi) frequency d/Sp of 
projections onto the Grassmanian manifold. Regarding (i) and (ii), we use Matlab's built-in Rungc- 
Kutta function ODE45 with adaptive step size control and set the absolute and relative tolerance 
to 10~ 6 and 10 -8 respectively. For (iii) the value of L^, at which the initial conditions for in 
(|3.32[) are specified, is selected so that \4>(x)\ < 10~ 14 for |x| > L M . The values for our 16 families are 



A n : 


Lao 


= 176, 


E on - 


Loo 


= 155, 


Aoff: 


Loo 


= 138, 


E oS : 


Loo 


= 135, 


Bon- 


Loo 


-75, 


F ■ 

A on ■ 


Loo 


= 270, 


-Boff: 


Loo 


= 103, 


F oS : 


Loo 


= 260, 


Con : 


Loo 


= 200, 


^on ■ 


Loo 


= 110, 


Cog- 


Loo 


= 140, 


GoB- 


Loo 


= 118, 


D on : 


Loo 


= 110, 


H on : 


Loo 


= 324, 


D oS : 


Loo 


= 125, 


H oS : 


Loo 


= 250. 



For (iv) we use the tolerance 10 10 on the L 2 norm of the residual of (|1.5|) . According to our 
tests |-E(0)| is insensitive to these parameters of (i-iv) if further refined (tolerances decreased and Loo 
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Fig. 3.2. Convergence of |E(0)| in dependence on the FEM discretization parameter h in the discretization of 
H1.5II and on the frequency d/8p of projections of the solution of H3.30I I onto the Grassmanian manifold. The selected 
SGS is one at T_ R3 0.48 in the family E a ff. In (a) we used d/Sp = 200 and in (b) h = 0.03. 

increased). We study next sensitivity to (v) and (vi). In both cases we observe a clear convergence 
of | £7(0 1 with refinement. Fig I3.2I shows dependence of |£/(0)| on h in the discretization of (| 1 . 5|) and 
on the projection frequency d/Sp for the case of an SGS at T_ m 0.48 in the family E q q with p = 3 
fixed. 

In our computations we use p = 3, h — 0.06, d/Sp = 50 as the default for the discretization of 
(|1.5p . We reduce h and/or increase d/Sp until \E(0)\ < 10 -8 . This threshold test cannot always be 
passed in the vicinity of folds where only 10 -6 is often achieved. In such cases it is checked that 
|£7(0)| < jq min Ae7 |i?(A)| for the selected contour 7 as an a-posteriori test. This criterion is satisfied 
in all the computations except when a zero of E lies in the vicinity of 7 for the particular SGS. 

Note also that with a finite difference (FD) discretization of (|1.5[) the values of |-E(0)| are much 
larger and neither convergence in dx nor in d/Sp has been typically observed, which has been tested 
with both 4th and 2nd order centered FD stencils. FD discretization is exptected to fail due to the 
profile 4> being only C 1 at the interface x = 0. 

3.5.2. Winding Number Computations. We present first winding number computations for 
the spectral problem as discussed in Sec. 13.31 In this section we do not attempt to determine the 
location of individual eigenvalues (if any) in the right half complex plane. 

The selection of the contour and its discretization is described in Sec. 13.31 We stress again that 
due to the small gap between the contour 7 and the imaginary axis we are able to capture only 
eigenvalues with Rc(A) > Rc(7). In all the computations in Fig. 13.41 and 13.51 Rcfo) = 0.005 was used. 
In Fig. 13.61 several examples with a smaller distance from the imaginary axis are also presented. 

Equation (|3.30[) was solved up to x = for both V + and V~ , i.e., up to the interface location. 
The Evans function was thus evaluated at = 0. After the evaluation of E along 7, the winding 
number n(E(^),§) was computed numerically using a simple Matlab script, see Appendix lAl 

For an SGS along the family A r Fig. 13.31 shows the curve E(^) extended symmetrically about 
the real axis using Lemma 13.21 The bold segments along E(-f) in part (d) correspond to segments 
of the contour 7 lying close to cr e ss(L) as marked in (c). Clearly, in the vicinity of <7 ess (L) the Evans 
function is more oscillatory than along the remaining parts of the contour. 

For the selected SGS families from Fig. 12.31 bifurcating from onsite GSs the results of winding 
number computations are in Fig. 13.41 Similarly, Fig. 13.51 corresponds to the SGS families from Fig. 
12.41 bifurcating from offsitc GSs. Full blue lines mark stable parts of the SGS families, where the 
winding number is zero, and the dashed red lines mark unstable parts, where the winding number 
is larger than 0. The winding number values are written next to the curve and correspond to the 
adjacent segment demarcated by the black circles. 
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Fig. 3.3. Evans function graph for a selected SGS. (a) The SGS family A jj, a selected SGS at T_ R3 2.063 
marked; (b) Profile of the SGS marked in (a); (c) the essential spectrum cr ess (L) (lie segments along Re(X) = 0) and 
the contour 7 (vertical dashed line) with segments close to cr ess (L) marked in bold red; (d) E(-y) for the SGS plotted in 
(b). Bold red segments mark the image of the marked segments along 7 in (c). 



Except for the cases A on , F on , A g and F a g each fold apparently marks a change in the value of 
the winding number. The reason why this does not happen for A on , F on , A Q g and F Q g is that in these 
cases all zeros of E have small real parts and when the zeros lie between the imaginary axis and 7, 
the resulting winding number is zero. Fig. I3.6I shows three of these cases for the contour 7 with 
Re(7) = 0.002, i.e. closer to the imaginary axis. The folds now more closely match winding number 
changes, which in these cases all correspond to stability changes. We show, however, in Sec. I3.5.3I 
computations suggesting that at least in the case of real eigenvalues bifurcations of eigenvalue from 
the imaginary axis happen near folds but not exactly at folds. In agreement with [41j we get that all 
off-site GSs (black dots in Fig. I3.5[) are unstable. For onsite GSs (black dots in I3.4|) we get 'stability' 
(Re(7) < 5 * 10~ 3 ) except for D on . Previously, see e.g [3TJ[T], onsite GSs, including discrete GSs, 
have been demonstrated to be stable via direct numerical simulations of their evolution in time. As, 
however, shown in 41, onsite GSs can, in fact, be unstable when an eigenvalue bifurcates from an 
edge of cr ess (L) away from the imaginary axis. This is caused by a resonance between the GS and an 
'internal mode'. From our selected cases D on is the only unstable one. Note that although the GS C on 
is plotted as stable in Fig. 13.41 it is shown in Fig. 6 of [U| that it is unstable due to an eigenvalue A 
with Rc(A) rj 3.5 * 10 -4 . This eigenvalue lies, however, well to the left of our contour. We reproduce 
eigenvalues near this point in the next section. 

3.5.3. Determining Eigenvalue Locations. Note that in general the method of the Evans 
function combined with the argument principle and a fixed contour cannot provide information on 
the strength of an instability since the precise location of eigenvalues is not known. When, however, 
the resulting winding number n is odd, the symmetry in Lemma 13.21 implies that at least one of the 
eigenvalues is real. Since in many of our examples n = 1, the corresponding instability is exponential. 
In these cases the eigenvalue can be located by evaluating E(X) throughout A € (0, ||r||oc||</>s||^o) and 



Surface Gap Solitons; Stability via the Numerical Evans Function 



25 



5 
4 

S 3 

UJ 2 
1 








- "stable" SGS (Re(*.)<0.005) 

- unstable SGS 



(2)1 <37 W 



0.5 

- 

50 


L 

-50 



5 r _ 

H 



10 



-1.5 



-0.5 



Fig. 3.4. Evans function winding number computations for selected SGS families bifurcating from onsite GSs. 
The vertical contour ■y with Re(y) = 0.005 was used. Full (blue) lines mark 'stable' SGSs, where n(_B(7),0) = and 
dashed (red) lines unstable SGSs, where n(E(-y),0) ^ 0. Changes in the values of n are marked by black circles and 
the value is written next to each segment of the curve. 



interpolating the graph to obtain roots of E. The graphs of the root A along the SGS family C on 
and along the family B Q g are plotted in Fig. I3.7I and l3.8l We see that in the case C on the eigenvalue 
does not bifurcate from precisely at the fold but shortly before the fold is reached, while for B Q g 
the bifurcation appears directly at the fold. 

For illustration in Fig. [3J]wc show the plot of E(X) along A G (0, HrHooH s ||So) for the SGS at 
the circle marker in the family C on in Fig. I3.7I 

In order to locate complex eigenvalues via the Evans function, there are, in principle, two al- 
ternatives: either a search algorithm employing the argument principle and contours enclosing each 
search region or a direct root finding approach based on a polynomial interpolation. We use the latter 
approach and select Miiller's method for finding complex zeros of complex functions [36] , Miiller's 
method is a generalization of secant method in that instead of a linear interpolation it uses a quadratic 
one and thus requires three points for an initial guess. We use the implementation of D.H. Cortes [8] 
and terminate Miiller's iteration if \E(\^)\ < 10~ 12 and |A^ _1 ) - A (fe) | < 1(T 6 for the fc-th iterate 
A' fe '. Selected results have also been checked with Matlab's standard routine fminsearch, which is 
based on a simplex search method. 

We restrict our attention to eigenvalues in or near the slit (0, 0.005) + iK omitted by the winding 
number computations in Sec. 13.5.21 and select two examples, the families C on and D g, for which we 
perform tracking of eigenvalues along the SGS families. For the C on gap soliton we select the starting 
guess triplet near 3.5 * 10~ 4 + 0.49 * i, which was reported in Fig. 6 [41] to be an eigenvalue. After 
convergence we track the eigenvalue throughout the whole SGS family. The results are in Fig. 13.101 
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Fig. 3.5. Evans function winding number computations for selected SGS families bifurcating from ofjsite GSs. 
The vertical contour 7 with Re(-y) = 0.005 was used. Full (blue) lines mark 'stable' SGSs, where n(E(-y),0) = and 
dashed (red) lines unstable SGSs, where n(E('y),0) ^ 0. Changes in the values of n are marked by black circles and 
the value is written next to each segment of the curve. 
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Fig. 3.6. Results of Evans function winding number computations for the SGS families F on , A g and F jj using 
the vertical contour^ with Re(-y) = 0.002. Apparently, a better match between the location of folds and stability changes 
is achieved than in the case Re(i) = 0.005, see Figs. \34\and[3H\ 



Varying the initial guess we have found two eigenvalues for the gap soliton: Ao ~ 7.8 * 10~ 4 + 0.498i 
and Ao « 0.5 * 10~ 4 + 0.4447 * i, satisfying our tolerances. The former eigenvalue is tracked in Fig. 
I3.10I (b) and (c) and the latter one in (d) and (e) . 

For the family D Q g and the starting guess triplet for the gap soliton near 0.002 + 0.544 * i the 
computed eigenvalues are plotted in Fig. I3.11I 

In both Fig. I3.10l and l3.11l the eigenvalue curves have a number of discontinuities. At these points 
Miiller's iteration switches to another eigenvalue. The computations show that the slit (0, 0.005) + iffi 
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Fig. 3.7. The real eigenvalue for the families C on and B g. Subplots (a) and (b) are for the family C on and (c) 
and (d) for B g. (a,c): SGS continuation curves; (b,d): the real eigenvalue along the family. Markers along the curves 
are for a better orientation. 
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Fig. 3.8. A zoomed-in view of the bifurcations of real eigenvalues from zero for the family Co 
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Fig. 3.9. The Evans function along (0, ||F|[oo ||<^s ||^o) for the SGS at the circle marker in the family C on in Fig. 




Fig. 3.10. Eigenvalue computation using Muller's method, (a) SGS family C on with several points labeled; (b) 
The real part of the eigenvalues Ao obtained by tracking the eigenvalues of the gap soliton (labeled '0'). The gap soliton 
eigenvalue was computed with the starting guess 7.8 * 10 — 4 + 0.498i. (c) The corresponding imaginary part, (d) the 
same as (b) except the starting guess for the gap soliton is 10~ 4 + 0.446 * i (at the curve ends in (d) the eigenvalue 
switches to that plotted in (b)); (e) the corresponding imaginary part. 
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Fig. 3.11. Eigenvalue computation using Muller's method, (a) SGS family D g with several points labeled; (b) 
The real part of the eigenvalues Ao obtained by tracking the eigenvalues of the gap soliton (labeled '0'). The gap soliton 
eigenvalue was computed with the starting guess 0.002 + 0.544 * i. (c) The corresponding imaginary part. 



missed by the winding number computations in many cases, indeed, contains additional eigenvalues. 
Numerically locating all these eigenvalues is possibly unfeasible with existing methods. 
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3.5.4. Direct Numerical PDE Evolution. For several selected SGSs we check the stability 
results by direct numerical simulations of the i-dcpendent NLS (jl.ip using the corresponding SGS 
profile </> plus a small perturbation as the initial condition. In detail, we take 



ip(x, 0) = (f>(x) + 0.01 max((f)(x)) sin' 



(^)sech 



where the second term represents an oscillatory perturbation with an expected large spectral support 
in <r(L) so that many modes of the linearization problem (|3.4|) are excited. This perturbation is 
used merely to accelerate possible instabilities since the numerical approximation alone introduces 
error which has a large spectral support. The numerical integration of was performed via a 

4th order (in dt) split-step method |51j where was rewritten as d t tjj — Aip + B{x, ip)ip with 

A = i9| and B(x,ip) = —iV(x) + iT(x)\ip(x)\ , so that splitting into the constant coefficient linear 
part d t ip = Aip and the nonlinear part d t ip = B(x,ip)ip was applied. The former part was solved in 
Fourier fc— space and the latter one in physical x— space exactly. In all simulations the computational 
box x <E [—100, 100 — dx] was used with 3000 grid points so that dx ~ 0.067. The time step was 
dt = 0.01. The two above mentioned conserved quantities of equation (|1.1|) . namely the total power 
IIV'C'j and the total energy E^ijj), are preserved with our numerical method with a typical 
relative error around 10 -5 and 10~ 4 respectively. 
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Fig. 3.12. Time evolution of 2 SGSs from the family Con (T+ = (a,) the SGS family with stability marked 
as in Fig. \3.j\ and 2 selected labeled points: [1] with T_ Ri 1.451 and [2] with T_ 1.472; (b) profiles of the 2 SGSs 
[1] and [2]; (c) and (d) evolution of the modulus \ij>\ for the 2 SGSs. Clearly, i/i^j evolves unstably and i/jp] stably in 
agreement with the stability diagram (a). 

Figs. 13. 1213.131 and [3441 show the evolution of one stable and one unstable SGS along the families 
C on ,G ff and F on , respectively. And Fig. 13.151 shows the evolution of the onsite GSs (r_ = T + = 1) 
-D on and E on . In all cases the (in)stability of the dynamics is in agreement with the results in Figs. 
13.41 13.51 and 13.61 In the case F on the time interval had to be chosen larger than in the other cases due 
to the weaker instability: for all the SGS in the F on family Re(A) < 0.005 while in the other families 
the unstable solutions have Rc(A) > 0.005 as follows from Figs. 13.41 13.51 and 13.61 

Based on a physical intuition from optics one might expect that SGSs with their center in the 
more defocusing half of the medium (i.e. the half with the smaller value of T) should be unstable due 
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Fig. 3.13. Time evolution of 2 SGSs from the family G g (T+ = — X). (a) the SGS family with stability marked as 
in Fig. 13.51 and 2 selected labeled points: [1] with T_ ss -1.764 and [2] with T- m —1.959; (b) profiles of the 2 SGSs 
[1] and [2]; (c) and (d) evolution of the modulus \ip\ for the 2 SGSs. Clearly, i/^jj evolves unstably and i/>[ 2 ] stably in 
agreement with the stability diagram (a). 




Fig. 3.14. Time evolution of 2 SGSs from the family F on (T_|_ = — X). (a) the SGS family with stability marked as 
in Fig. \3.4\ and 2 selected labeled points: [1] with T_ -1.643 and [2] with T_ Si —2.229; (b) profiles of the 2 SGSs 
[1] and [2]; (c) and (d) evolution of the modulus \vb\ for the 2 SGSs. Clearly, ip^ evolves unstably and i/)[ 2 ] stably in 
agreement with the stability diagram (a). 



to the tendency of light to converge to the more focusing regions. Mathematically this expectation 
also makes sense since the total energy (11.61) is decreased via a shift toward the more focusing region. 
This is, however, only a heuristic view due to the lack of shift invariance in the system. In reality 




Fig. 3.15. Time evolution of the GSs D on and E on (T+ = T_ = 1). In Fig. \3.4\ Evans function computations 
show these to be unstable and stable respectively. The t— evolution agrees with these stability results. 



nothing, in general, prevents the possibility of existence of a local extremum of the energy with the 
extremizer localized in the less focusing region and being stable. As we can see, for instance, in the 
SGS families C on and F a g by inspecting Figs. 12.81 and 13.41 there are parts of the family where this 
intuition indeed fails. The segments of the families marked in Fig. 13.161 correspond to stable SGSs 
centered in the less focusing half. Note that for the segment of the F g family exiting the plot window 
on the right the part that satisfies T_ > is a family of "stable" solutions (Re(A) < 0.005) at a 
focusing/dcfocusing interface and centered in the dcfocusing medium! Similar segments of families of 
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Fig. 3.16. Left and middle: "stable" segments (Re(X) < 0.005,) of the SGS families C on (where r+ = I) and 
F a ff (where T+ = —1), which consist of solutions centered in the less focusing part of the domain are marked in bold. 
On the right: profile of the "stable" SGS (1) at V— ~ 1.41 marked along the F a ff family. This SGS is at a truly 
focusing /defocusing interface and concentrated mostly in the defocusing half. 

SGSs centered in the less focusing half can be found, for instance, in the cases A on , F on , A s, and C g. 

4. Conclusions. Surface gap solitons (SGSs) of the cubic ID periodic Schrodingcr (PNLS) / 
Gross Pitaevskii (GP) equation with a nonlinearity interface have been studied. As the computa- 
tions show, bright SGSs exist even in the case of a focusing/dcfocusing interface. In [133 SGSs were 
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constructed via a numerical parameter continuation from gap solitons (GSs) and the families were 
followed up to the first fold. Here we have extended these families using the numerical arclcngth 
continuation. As shown, the SGS families often contain GSs centered at different cxtrema of the 
linear periodic potential or several non-identical GSs at the same extremum. For families of spatially 
broad SGSs the number of folds appears to equal the number of shifts of the SGS center of mass from 
one extremum of the periodic potential to the next. Asymptotic methods reveal that for SGSs bifur- 
cating from GSs centered arbitrarily far from from the interface the allowed size of the jump in the 
nonlincarity coefficient decays as the propagation constant /frequency approaches the corresponding 
bifurcation edge of the spectral gap. Numerical arclength continuation results confirm this statement. 

To inspect linear stability of SGSs we have developed a numerical Evans function method. Eigen- 
values with the real part to the right of a small threshold have been detected via the complex argument 
principle applied to the Evans function with a vertical contour near the imaginary axis. This is the 
first time, to our knowledge, that the Evans function method combined with the argument principle 
has been applied to a problem with periodic coefficients. Standard difficulties in evaluating the Evans 
function in a numerically stable way have been reported and explained in detail in the context of 
the problem at hand. These include firstly stiffness of the ODE linearization problem on each of the 
unstable and stable manifolds, which need to be approximated for the Evans function evaluation. This 
is overcome via a reformulation in exterior algebra. Secondly, a numerical approximation of solutions 
with a fast exponential growth is necessary, which is performed under a change of variables. Next, 
accurate winding number computations require an adaptive discretization of the contour. Finally, the 
accuracy of the Evans function is shown to depend on the level of preservation of the weak Grassma- 
nian invariant in the evolution of the two manifolds and on the discretization error in the solution of 
the SGS profile. 

The results reveal the existence of both unstable and stable SGSs, where stability does not exclude 
the presence of eigenvalues with a positive real part below the threshold given by the distance of our 
contour to the imaginary axis. The results on stability of GSs are in agreement with the previous 
results in [41j . Interestingly, stability is possible also for SGSs centered in the less focusing medium 
of the two and even when the SGS is concentrated in the defocusing half with a focusing-defocusing 
interface. Direct numerical simulations of the PNLS/GP equation were performed for several SGS 
examples and confirm the (in)stability statements based on Evans function computations. 

An apparent shortcoming of the argument principle approach is the lack of information on eigen- 
values between the imaginary axis and the vertical contour. Note that the imaginary axis cannot be 
chosen as the contour even if the Evans function E is analytically extended through the essential spec- 
trum, see [231 [231 [35] for an analytic extension in constant coefficient systems, due to the occurrence 
of zeros of E on the imaginary axis. In order to gain insight on the region between the imaginary axis 
and the vertical contour, we have tracked several complex eigenvalues there via the use of Miiller's 
method. The authors plan to address numerically eigenvalue bifurcations from the essential spectrum 
by using an analytic extension of the Evans function in a future paper. 

All the techniques presented in this paper can be easily adapted also for problems with linear 
interfaces, e.g. when the linear potential has the form V(x) = V+(x) X[o,oo)( x ) + V-i, x ) X(-oo,o){ x ) 
with periodic V± . In fact, only periodicity of V for large enough |a;| is a requirement for the techniques. 
The stability method can thus be directly applied, for instance, to the linear interface SGSs in [2"51 
134] . In addition the form of the nonlincarity coefficient can be arbitrary. A practical advantage 
of the numerical Evans function method is the possibility of a straightforward parallelization of the 
computations by dividing the contour into segments. 

Computation of SGSs via the arclength continuation generalizes directly to more spatial dimen- 
sions. The application of the Evans function is, however, generally unfeasible as the stable and 
unstable manifolds are typically infinite dimensional in more than one spatial dimension. 
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Appendix A. Winding Number Code. For readers' convenience and to motivate the use of 
winding number based computations of point spectra by others, we provide here our Matlab code for 
the computation of the winding number of a closed curve in the complex plane with respect to the 
point 0. A remark on the choice of a starting point is included at the beginning of the code. 



function wind = wincLnr (E) 



% winding number computation 

% input: complex array E of the size [N x 1] representing a closed curve 

% in the complex plane ( E(l) = E(N) ) 

% output : winding number wind of the curve with respect to the point 



% beware: choice of ind_st 



plot s_on = l ; 
if (plots.on) 

figured); plot(E) 

xlabel ( 'Re (E) ' ) ; ylabel ( ' Im (E) ' ) 

end 

%select a starting point for the winding number computation 

% (should not be at a point where the curve E is tangent to a ray from the origin 
% or where the phase is \pm \pi) 
ind_st=3 ; 

%for the Evans function example with an unbounded lambda— contour where the 

%first m and the last m entries of E are a line connecting the last computed value 

%of E along the lambda— contour and the asymptotic value 1 of E for 

% | lambda|— > oo a suitable choice of ind_st is 1 < ind_st < m because the complex 

%phase of E changes monotonically and the phase value is not +/— pi for Kind_st<m 



%initiate the value of the phase and the winding number 
phase.O = atan2 ( imag (E ( ind.st )), real (E ( ind.st ))) ; 
phase = atan2 ( imag (E ( ind.st + 1 ) ) , real (E ( ind.st + 1 ) ) ) ; 
if (ind.st==l) 

phase.old = atan2 (imag (E (end— 1) ), real (E (end— 1) )) ; 

else 

phase.old = atan2 ( imag (E ( ind.st — 1 )), real (E ( ind.st — 1 ))) ; 

end 

if (phase>phase_0 && phase_0>phase_old && abs (phase— phase.old) <pi ) 
wind = 1; 

elseif (phase<phase_0 && phase_0<phase_old && abs (phase— phase.oldXpi ) 
wind = —1; 

end 

%run along the E— curve tracking the phase and increasing or decreasing 
%the winding number counter 
check_monot = 0; 
phase.old = phase.O; 

for k = [ind.st + 1 : length (E) 2:ind.st] 
phase=atan2 (imag (E (k) ) , real (E (k) ) ) ; 
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if (check_monot==l ) %check if E is tangent to the ray with the angle phase_0 

if (phase > phase.old && phase.olcLold < phase.old && abs (phase— phase_old)<pi) 
wind = wind + 1 ; 

elseif (phase < phase.old && phase_old_old > phase.old && abs (phase— phase_old)<pi) 
wind = wind — 1; 

end 

check_monot = 0; 

end 

%check if the value phase_0 has been passed 

%(i.e. if E has crossed the ray with the angle phase_0) 

if (phase>phase_0 && phase_old<phase_0 && abs (phase— phase_old)<pi ) 

%counter clock— wise crossing 

wind = wind + 1; 

elseif (phase<phase_0 && phase_old>phase_0 && abs (phase— phase_old)<pi) 
%clock— wise crossing 
wind = wind — 1; 

elseif (phase==phase_0 ) %need to check if phase has an extremum here 
%(i.e. if E is tangent to the ray with the angle phase_0) 
phase_old_old = phase.old; 
check_monot = 1; 

end 

phase.old = phase; 

end 
end 
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